Evolutionary dynamics, evolutionary forces, and robustness: A nonequilibrium statistical mechanics perspective

Significance Evolution through natural selection is an overwhelmingly complex process, and it is not surprising that theoretical approaches are strongly simplifying it. For instance, population genetics considers mainly dynamics of gene allele frequencies. Here, we develop a complementary approach to evolutionary dynamics based on three elements—organism reproduction, variations, and selection—that are essential for any evolutionary theory. By considering such general dynamics as a stochastic thermodynamic process, we clarify the nature and action of the evolutionary forces. We show that some of the forces cannot be described solely in terms of fitness landscapes. We also find that one force contribution can make organism reproduction insensitive (robust) to variations.

A conventional view of evolutionary dynamics is based on three essential elements (1): 1) organism reproduction with imperfect heredity; 2) variations, including mutations, which are typically introduced by the reproduction process; and 3) selection, which acts within a population and allows some variant species to survive and reproduce, while eliminating others. When considering variations, a sizeable fraction of evolutionary biology is focused on genetic and epigenetic variations. However, variations upon which selection acts are occurring on multiple levels, and involve many entities, traits, and behaviors that are usually encapsulated by a rather imprecise concept of phenotype (2,3). Regarding selection, many phenotypic aspects contribute to long-term survival and reproduction. Two instances are the interactions between the organisms (e.g., sexual reproduction, predation, competition and cooperation, and social organization) and the interactions with biotic and abiotic environmental factors (e.g., the presence of other species or the inorganic composition of a certain habitat), whose changes span multiple spatiotemporal scales.
The main role of (necessarily) simplistic mathematical models of evolution is to analyze the possible outcomes of evolution and to explore the assumptions that generate these outcomes. In this way, one hopes to clarify some essential concepts used by evolutionary narratives. In the present work, we formulate a simple model that incorporates the three essential elements described above. Consequently, we call the evolutionary dynamics described by this model the reproduction-variation-selection (RVS) dynamics. In formulating the model, we seek both simplicity and generality. For the sake of generality, we specify neither the particular nature of the hereditary variables nor that of the related variations: They can be genetic, epigenetic, or phenotypic. For the sake of simplicity, we assume large population sizes and the presence of a constant environment. By building a theoretical framework inspired by nonequilibrium statistical mechanics (4-7), we can study this dynamics in its generality. In fact, we can clearly define the notion of evolutionary force and explicate its relation to reproduction, variations, and selection. Because of the simplicity of the model, we are able to get analytical expressions for different force terms, and to perform explicit analyses of the role played by them. In particular, we uncover an evolutionary force within the RVS dynamics that can engender robustness of reproduction to variations, without any explicit selection for this trait. By adding restraining assumptions to the model, we can also make simple predictions about the behavior of population robustness during the RVS dynamics. Lastly, we compare our predictions with the results of laboratory experiments on populations of evolving viruses.

Model
Generic Evolutionary Dynamics. Generic evolutionary dynamics can be described as follows. We imagine a population of N organisms, each of which is characterized by some hereditary variables γ. These variables can describe genetic and epigenetic factors, collective phenotypic features, etc. We do not specify the precise biological nature of γ, but Significance Evolution through natural selection is an overwhelmingly complex process, and it is not surprising that theoretical approaches are strongly simplifying it. For instance, population genetics considers mainly dynamics of gene allele frequencies.
Here, we develop a complementary approach to evolutionary dynamics based on three elements-organism reproduction, variations, and selection-that are essential for any evolutionary theory. By considering such general dynamics as a stochastic thermodynamic process, we clarify the nature and action of the evolutionary forces. We show that some of the forces cannot be described solely in terms of fitness landscapes. We also find that one force contribution can make organism reproduction insensitive (robust) to variations.
we assume that they determine the reproductive success of each organism. Mathematically, γ can be represented as multidimensional continuous or discrete variables, or a mixture of the two; here, for simplicity's sake, we focus on the discrete case (but see SI Appendix, section II for the continuous case). We refer to an organism characterized by the variables γ as a being of type γ.
The typical number of offspring that each type engenders in one generation is the reproduction rate f γ . When the environment is assumed constant, selection favors organisms reproducing faster. During one generation, random changes of the hereditary variables create the diversity upon which natural selection can act. We generically refer to these changes as variations. For the particular case of genetic variations, they may involve a single organism, or pairs of organisms, through recombination. As stated in the Introduction, variations can also be epigenetic and/or phenotypic.
Over many generations, the evolutionary dynamics can be described as a discrete-time stochastic Markov jump process. The population composition is denoted by the vector n = (n γ ): each of its entries represents the number of γ-type present in the population. The probability of observing a certain n at generation τ , p n (τ ), is described by a Chapman-Kolmogorov equation (see Box 1 for a descriptive summary of the main formulae), p n (τ + 1) = n W n←n p n (τ ). [1] The entries of the transition kernel W n←n determine the probability that the population transitions from n to n over one generation. Notice that we assume nonoverlapping generations. The precise form of W n←n depends on the details of the stochastic evolutionary process. In general, f γ , and thus W n←n , changes with the environment. Before committing ourselves to a particular W n←n , let us introduce a quantity that we call evolutionary directionality. It is defined as the log ratio of forward and backward transition probability, Evolutionary directionality will play a crucial role in our discussion. For a pair of transitions n n such that W n←n W n ←n , neither direction is probabilistically favored. This situation corresponds to negligible directionality, F n←n 0. In contrast, when n n are such that W n←n W n ←n (respectively, W n←n W n ←n ), the population is more likely to evolve in the direction n ← n (respectively, n ← n). This situation corresponds to nonzero directionality, F n←n > 0 (respectively, F n←n < 0).
In general, the evolutionary directionality F n←n can be written as a sum of two contributions, the first of which can be expressed as a potential difference, while the second cannot, The precise form of ψ n and ζ nn depends on how selection and variations affect the population, expressed in the model through W n←n . We refer to ζ nn and ψ n − ψ n as evolutionary force contributions, and they are important for the following reason. If ζ nn is positive for n ← n , then ζ nn positively contributes to the directionality F n←n , and we can say that this force contribution +const., [11] where μ γ and ω γ are given in Eqs. 4 and 9.
Restricted RVS Dynamics-Additional assumptions: -Strong selection (SI Appendix); -Type-independent and small variation Prob. μ γ μ. Expected Growth potential where F n := γ f γ n γ and Ω n is given in Eq. 14.
favors such transition. Analogously, if ψ n increases along n ← n , then ψ n − ψ n is positive, and such a transition is favored by the force contribution originating from ψ n . In analogy to similar quantities encountered in nonequilibrium statistical mechanics and thermodynamics (see, e.g., ref. 8), we refer to the two contributions in Eq. 3 respectively as conservative and nonconservative evolutionary forces. Note that, when the latter force vanishes, ζ nn = 0 for all n and n , one can show thatp n ∝ exp {ψ n } is the stationary distribution of the dynamics described by [1] (see, e.g., ref. 9). This situation corresponds to a gradient stochastic dynamics in which the landscape identified by ψ n determines the conservative (or gradient) force ψ n − ψ n .
RVS Dynamics. A model of evolutionary dynamics is specified by a particular choice of the kernel W n←n . In our RVS dynamics, this kernel is defined through the probability that a variation γ ← γ occurs over a single generation, π γ←γ (where π γ←γ ≥ 0 and γ π γ←γ = 1 for all γ ). We consider reversible variations, so that π γ←γ > 0 ⇔ π γ ←γ > 0. This assumption does not preclude that variations in some of the directions may be much more likely than in the opposite ones, that is, π γ←γ π γ ←γ for some γ and γ . Crucially, π γ←γ must be regarded as a function of γ . Indeed, π γ←γ is also subject to evolution: Different types may evolve to vary in different ways. In general, many detailed mechanisms can contribute to variations of a type, all encapsulated in the probability π γ←γ . In addition, for variations involving pairs of organisms, the corresponding contribution to π γ←γ depends on the population composition n. For the sake of simplicity, here we do not distinguish among different mechanisms of variation (see SI Appendix, section IE for such details). We consider only the overall probability that γ has varied, [4] In our model, the creation of a new pool of variants and the selection from this pool take place in each generation, as in Fig. 1. We assume the overall size of the population to be large, N = γ n γ 1, so that its fluctuations can be neglected and its value regarded as a constant SI Appendix, section IB. This might reflect the fact that environmental resources are abundant but limiting. For simplicity's sake, we also assume that the organisms do not directly interact, so that the reproduction rate f γ does not depend on the population composition n. The case of interacting populations with small and fluctuating population sizes is analyzed in SI Appendix, section IA.
Under these assumptions, the transition kernel can be expressed in a simple form (details in SI Appendix, sections IA-IC), [5] Note the alternation of variations, reproduction, and stochastic selection, as in Fig. 1. The term in square brackets describes 1) the probability that variations affect the organisms, γ π γ←γ n γ , and 2) organismal reproduction, f γ . Then, the multinomial prod- nγ ensures the random selection of N organisms.
Fast-growing types-types with larger numbers of descendants, f γ -are more likely to be selected, yet selection fluctuations (usually referred to as genetic drift in the literature of evolutionary dynamics) are accounted for. The model [5] can be regarded as a generalized version of the Wright-Fisher model in which 1) Schematic representation ofthe evolutionary dynamics as modeled by the transition probabilities in Eq. 5. The symbols •, , and (circle, star, and square) represent different types of γ present at generation τ − 1, in a population of N = 7 organisms. As variants are generated-dashed linesa new type appears, . Continuous lines denote no variation. The variation probabilities π γ→γ as well as the probability that no variation occurs, 1 − μγ , are also reported. Selection-finely dashed lines on the right-finally determines which types and in what amount make it to the next generation. The reproduction rate fγ determines the chances of being selected.
hereditary variables, as well as their variations, are generic and 2) organisms are subject to variations, and then selected preferentially by their reproduction rate (cf. ref. 10, equation 3.68).

Results
Evolutionary Forces. For the model described by the kernel [5], the conservative force potential and the nonconservative force contribution are where φ n := N ln γγ f γ π γ←γ n γ .

[8]
As discussed in SI Appendix, section IC, these expressions follow from a simple algebraic procedure using [5]. This general procedure is justified by the fact that we do not specify the nature of γ, nor how variations affect γ.
We can now analyze the expressions for ψ n and ζ nn and find conditions for which ζ nn is positive, and for which ψ n increases. As mentioned before, such conditions indicate compositions n toward which the population is more likely to evolve.
The nonconservative force ζ nn accounts for the effect of variations and selection fluctuations (aka genetic drift). This fact is highlighted by the presence of the variation probabilities π n←n and the multinomial product in Eq. 7. Although, in general, one cannot determine easily the sign of ζ nn , two limiting cases are useful to gain some insight into the effect of this nonconservative force (details in SI Appendix, sections III and ID). We first consider the limit of small variation probabilities [4], μ γ μ → 0. Imagine a population composition transition n ← n in which a particular type γ disappears: n γ = 0 ← n γ = 0. This event is irreversible, since the nonconservative force diverges logarithmically (ζ nn = −ζ n n ≈ − ln μ → ∞, for μ → 0), thus preventing γ reappearance. Generalizing this insight to small but finite μ γ , we can argue that ζ typically causes a loss of diversity by favoring the selection of the most abundant types. Second, consider the idealized case when the variation probability solely depends on the variant type, that is, π γ←γ = π γ for all γ and γ . The probability π γ quantifies the likelihood that γ is engendered by variations, irrespective of the type from which it originates. In this case, the nonconservative force actually becomes conservative, ζ nn = γ n γ ln π γ − γ n γ ln π γ . It generically favors n with higher γ n γ ln π γ ; that is, it favors population compositions with a higher cumulative log likelihood of being engendered.
The potential ψ n that appears in the conservative force is composed of three terms (Eq. 6). The first term is associated with selection, since it is simply the cumulative logarithmic growth of the population. The corresponding evolutionary force, γ (n γ − n γ ) ln f γ , drives the population toward compositions with a higher cumulative log reproduction rate.
The second term is an entropic contribution, and it arises due to the assumed indistinguishability of organisms of the same type. The corresponding force, γ ln n γ ! − γ ln n γ !, favors heterogeneous population compositions. As a limiting case, consider again a transition in which a particular type γ disappears: n γ = 0 ← n γ = 0. Such an event is disfavored by this entropic force, which scales as the logarithm of the population size (∼ − ln N < 0). Hence, for constant μ and μN > 1, the entropic force counteracts the loss of diversity generated by the nonconservative force ζ (details in SI Appendix, section III and Eq. S51).
Finally, the last term in Eqs. 6 and 8 can be associated with the combined effect of selection and variations. We name this term expected growth potential of the population n, since it quantifies the expected growth of n at the next generation. It is larger for populations that are likely to generate variants with a high reproduction rate. This kind of population is thus favored by the force term Since φ n depends on n in a nonlinear fashion (see the logarithm), φ n can be viewed as describing an effective interaction acting among the types present in the population. This interaction is engendered by the combined action of selection and variations. All main mathematical expressions are summarized in Box 1, whereas the nature and effect of each evolutionary dynamical force ( [6] and [7]) are illustrated in Box 2 using two simple toy models.

Robustness.
Robustness is an important quantity that measures the extent of insensitivity of an organism, or its particular phenotypic features, to extrinsic or intrinsic variations. For instance, the so-called mutational robustness measures how much a given phenotypic feature changes upon genetic mutations (in a given environment). The expected growth potential φ n can be explicitly related to robustness. To see this, let us first introduce ω γ , a measure of the sensitivity (i.e., the inverse of robustness) of the reproduction rate to type variations, Types γ with sufficiently large reproduction rate f γ have positive sensitivity, ω γ > 0; see also 1) imply that the relative decrease of the reproduction rate, (f γ − f γ )/f γ , averaged over all possible variations, π γ ←γ /μ γ , is large. Hence, types with large ω γ have large sensitivity, or small robustness to variations. In contrast, ω γ 0 characterizes those types whose reproduction is insensitive-and hence robust-to variations.
The expected growth potential φ n can be now expressed as a decreasing function of ω γ , that is, the lower the sensitivity of the organisms of the population, the higher φ n ; see SI Appendix, section IE.
To gain further understanding of Eq. 10, let us consider populationsn composed of types with relatively high and roughly equal reproduction rates: f γ f high , for any γ such thatn γ > 0. This case could describe situations in which selection is so strong that only the fastest-reproducing types survive. Eq. 10 can be approximated as φn − γ μ γ ω γnγ + const. [11] (SI Appendix, section IE). The right-hand side can be interpreted as the cumulative variation probability times sensitivity of the populationn. Hence, there are only two ways in which the potential term φn can increase: Either the types of the population can be subject to less variations, that is, μ γ decreases, or they can be more robust to variations, that is, their sensitivity ω γ decreases. Eqs. 10 and 11 have thus a general, important significance: They show how the evolutionary force generated by φ n − φ n favors populations whose types exhibit either low variation rates, μ γ 1, or high robustness to variations, ω γ 1. Crucially, the emergence of these features is a property of generic RVS dynamics and is not restricted to either specific hereditary variables or specific variation mechanisms.
Restricted RVS Dynamics and Fluctuation Relation. Eqs. 3, 6, and 7 specify the forces driving the evolution of RVS dynamics and how these forces shape its evolutionary outcome. However, the stochastic dynamics described by Eqs. 1 and 5 cannot be solved exactly. Hence, in general, it cannot be rigorously assessed to which extent each force contributes to the evolutionary outcome. To partially overcome this limitation, we will now discuss a simple yet biologically relevant case, for which we do have an approximate analytical solution.
First, we neglect the nonconservative force ζ nn [7]. Heuristically, this can be justified when the variation probabilities μ γ are not too small and selection is strong compared to selection fluctuations-types are highly discriminated by their reproduction rate; see SI Appendix, section IV. These conditions guarantee that ζ nn remains finite and that the effect of the forces that involve selection predominates over ζ nn (we recall that ζ nn may diverge for small μ and does not account for reproduction and selection). When neglecting ζ nn , the dynamics becomes approximately conservative, and p n (τ ) converges to a Boltzmann-like probability distribution-provided that W n←n [5] is ergodic- where Z := n exp ψ n is a normalization factor. For long times, population compositions with higher ψ n are exponentially more likely to be observed. Second, we limit our analysis to a regime in which the variation probability is small and similar for all types: μ γ := γ =γ π γ ←γ μ, for all γ. We can hence approximate the conservative potential [6] as (details in SI Appendix, section IV) where F n := γ f γ n γ is the cumulative growth, and with F πn := γ,γ ( =γ) f γ π γ ←γ n γ /μ being the cumulative growth upon variations. In analogy with the sensitivity of the type ω γ , Eq. 9, Ω n can be viewed as the sensitivity of the whole population. Ω n is also akin to so-called variation loads, where F n is replaced by some maximal or reference population reproduction rate (see, e.g., refs. 11 and 12).
In the asymptotic equilibrium regime [12], the population behaves similarly to a thermodynamic system in equilibrium with the environment. The population sensitivity Ω n can be regarded as an energy function, and μN can be regarded as an inverse temperature. Since ψ n depends, in a nonlinear way, on n, the detailed properties of the distribution [12] remain difficult to assess. However, a simple calculation inspired by equilibrium statistical mechanics (see ref. 13 and SI Appendix, section IV) allows us to gain insights about how p n (∞) affects the population sensitivity Ω n . Indeed, where Ω n = n p n (∞)Ω n denotes the average sensitivity over the equilibrium distribution. This relation binds the fluctuations of the sensitivity (right-hand side) to the changes of its average (left-hand side). In the context of thermodynamics, d Ω n /dμ are called susceptibilities, as they quantify how a certain macroscopic observable, here Ω n , is susceptible to the changes of a certain parameter, here μ. Eq. 15 tells us that the fluctuations of Ω n are higher when Ω n is more susceptible to changes of μ.
There are three important implications of the fluctuation relation [15]. First, the average population sensitivity Ω n decreases as a function of the variation probability, d Ω n /dμ 0. Since the inverse of Ω n measures the population robustness to variations, this result shows that robustness increases with μ. Second, the fluctuations of sensitivity vanish for N → ∞. Third, assuming that Ω n saturates at high values μ reaching a plateau, the fluctuations Ω 2 n − Ω n 2 are expected to decrease with μ. The fluctuation relation [15] is the main result obtained with additional assumptions imposed on the general RVS model. Although this relation has been derived with additional assumptions and is only approximate, it seems consistent with published results of laboratory experiments performed on viral populations (see Discussion and Conclusions and SI Appendix, section VIII). Comparisons with numerical simulations are discussed in Box 2 and in SI Appendix, section V.

Discussion and Conclusions
From a physicist's point of view, evolution is an example of a nonequilibrium stochastic dynamical system (14)(15)(16)(17)(18)(19). One of the main reasons why it is very different from dynamics studied in physics or chemical physics is that the reproduction of organisms is tightly connected with relatively precise, long-time inheritance [in contrast to self-reproduction of simple molecular systems, such as amphiphilic micelles (20)]. This allows the existence of stable variants that can better survive and reproduce than others, and thus allows organisms to evolve. Since the so-called evolutionary Modern Synthesis (21,22) was deeply anchored in population genetics (23), mathematical models of evolutionary dynamics have mostly dealt with genetic variations. Evolutionary Forces Clarified. In this work, we have tried to move away from this dominant scheme of modeling of evolutionary processes. Our RVS dynamics model is closely related to stochastic models of nonequilibrium thermodynamic processes (4)(5)(6)(7). There are no assumptions about the nature of either the hereditary variables or the variations: They can be genetic, epigenetic, or phenotypic. The attractive aspect of our simple model is that it clarifies the notion of "evolutionary forces." These forces emerge from stochastic variations and selection. They can be divided into two classes: "conservative" ones, that can be expressed in terms of potentials, and "nonconservative" ones, that cannot, Eq. 3 (8,24). Nonconservative forces, well-studied in many physical systems, have not been explored until recently in evolutionary biology. Their effect is not intuitive, since a usual "landscape metaphor" cannot be easily applied to them (25). This is not the case for the conservative forces, which can be easily intuited. In our model, we find three separate terms which contribute to the conservative force Eq. 6. The first term simply describes selection for fastest-reproducing variants, while the second is purely entropic and induces population diversification. The remaining third conservative force contribution can be related to reproduction robustness of individual types, Eq. 10. The existence of this "robustness-generating" evolutionary force is, in fact, a somewhat surprising result of our analysis.
Our dynamical analysis of evolutionary forces also clarifies an academic controversy about the interpretation of these forces (26). It was indeed debated whether selection, variations, and selection fluctuations should be interpreted as Newtonian forces (namely, causes of changes), or as statistical pseudoforces (i.e., stochastic events which affect the population). From our analysis, it becomes clear that these three elements should be interpreted

Box 2. Illustrative Toy Models
We consider here two examples whose purpose is to illustrate the expression of the evolutionary forces, the enaptation of robustness, and the fluctuation relation.
Two-types Model: Nature of Evolutionary Forces. Consider an RVS dynamics with three organisms, N = 3, and two types, γ 1 and γ 2 . Both types vary with the same probability, π γ1←γ2 = π γ2←γ1 = μ ≤ 1/2, and-without loss of generality-f γ1 > f γ2 . The network of type variations and that of transitions in population space are In the latter space, the reproduction rate and entropic forces act respectively as follows (details in SI Appendix, section VI): (2, 1) where the number of arrows is proportional to the strength of the force. The reproduction rate force favors the population with the highest proportion of fast-reproducing types, (n 1 , n 2 ) = (3, 0), while the entropic force favors a diverse population composition, (n 1 , n 2 ) ∈ {(1, 2), (2, 1)}. Note that, since both of these forces are conservative, the sum of their values along any cycle vanishes. In contrast, the nonconservative force favors homogeneous compositions, and the sum of its values along cycles does not vanish, in general, This force is stronger on the diagonal transitions (double arrows) than on the other ones. This can be explained by the fact that the loss of types ("loss of diversity") is higher: Two individuals of a certain type are lost rather than one (vertical transitions, single arrow) or none (horizontal transitions, no arrow). The right panel depicts how this force acts along cycles. See SI Appendix, section VI for mathematical details and a discussion about the force arising from the expected growth potential.

Square Grid Model: Enaptation of Robustness.
To illustrate the effect created by the force arising from the expected growth potential, Eqs. 8 and 11, consider the model depicted below, in panel a. Each site of a square 12 × 12 grid represents a distinct type, and variations consisting of changes of one type into another correspond to transitions between nearest neighbors. We assume the variation coefficient to be constant and equal toμ. Islands of fast-reproducing types (big orange sites) are surrounded by a sea of slowly reproducing ones (small gray sites), and each island differs by connectivity, namely, the number of adjacent slowly growing types. Among types with the same reproduction rate, the expected growth potential favors those surrounded by fast-reproducing types, namely, those belonging to the interior of the top left island. From numerical simulations, the steady-state type distribution (right heat map plot, logarithmic scale) clearly shows a higher probability of finding types surrounded by equally fit types. Indeed, 1) the probability of finding an isolated fast-reproducing type (isolated orange sites) is < 10 −10 , and 2) the probability of finding types in the upper square is roughly four orders of magnitude higher than that of finding types in the lower one-dimensional strip. This model illustrates how robustness, that is, insensitivity to variations, is enapted by the evolutionary dynamics.

Square Grid Model: Fluctuation Relation.
To illustrate our fluctuation relation and its implications, we use additional numerical simulations of the grid model. The plot below depicts the average population sensitivity, Ω n vs. the variation probability log 10 μ. The vertical bars represent 1 SD, Ω 2 n − Ω n 2 . This plot demonstrates the good qualitative agreement between the fluctuation relation and numerical results. As predicted by [15], 1) Ω n decreases as μ increases, and, 2) as Ω n approaches the plateau, fluctuations decrease as well. We refer to SI Appendix, section V for a more detailed discussion. as pseudoforces, but, nonetheless, they combine into Newtonian forces, whose expression appears in Eqs. 3, 6, and 7.

Robustness to Variations as an Enaptation.
Robustness is a concept that has attracted a lot of attention in biology (27)(28)(29)(30)(31)(32). 1) Many phenotypic traits of organisms and their developmental pathways have been shown to exhibit robustness to genetic mutations (33), intracellular component concentrations (34), external perturbations (35), etc. There have been also observations of robustness in the behavior of whole multispecies ecological systems (36,37). Finally, life on Earth has been, until now, robust enough to sustain life for billions of years despite strong internal variations and external perturbations. 2) Many studies of model organisms were performed in laboratories to uncover the molecular and cellular mechanisms which underlie robustness (38). Generic mechanisms-such as redundancy, modularity, and feedback control-have been described in many systems (39)(40)(41)(42)(43).
3) It is difficult to demonstrate, but it is widely believed, that robustness facilitates better survival and reproduction in naturally occurring biological systems (44). 4) It is possible to show, using computer simulations and simple mathematical models, that, under some important assumptions (e.g., strong selection, and high mutation rates in sequential genomes), robust variants may outcompete even faster-growing variants (45)(46)(47)(48).
The general mechanism behind this last phenomenon-which has been referred to as "survival of the flattest"-is elucidated by our analysis. In fact, the existence of a "robustness-generating" evolutionary force points toward a potentially novel mechanism associated with evolution-a mechanism that falls under the concept of neither adaptation nor exaptation (49). * A phenotypic feature, growth robustness, does not emerge in our model as a consequence of selection for robustness because of its fitness value, nor is it coopted from a previously evolved trait for its new function. Rather, the "robustness-generating" evolutionary force is intrinsic to the evolutionary dynamics itself. This force contribution expresses the tendency of the selected variant to be "surrounded" (in the type space) by other similarly fit variants. We propose to call this class of emergent phenomena "enaptations" (from Latin en-, in, into, + aptare, to fit), since robustness emerges here as a consequence of intrinsic aspects of evolutionary dynamics per se.
We would like to stress again that the emergence of robustness through enaptation does not preclude the existence of other mechanisms which engender this trait (see, e.g., ref. 29).
Quasispecies, Viruses, and a Fluctuation Relation. It might have not escaped the attention of the reader that our model has some similitude with previous models of evolutionary dynamics like the Wright-Fisher or the quasispecies models (10,(51)(52)(53); see SI Appendix, section VII. In contrast to these models, the nature of γ, as well as their variations, are generic. Quasispecies models belong to the class of our restricted RVS dynamics, since they typically consider genetic types with sequential genomes of finite length, and limit variations to point mutations, whose overall probability is the same for all types. More interestingly, the same assumptions seem to be approximately satisfied for many laboratory experiments studying evolution of virus populations * Exaptation (from Latin ex-, out of, + aptāre, to fit) refers to the recruitment of a biological trait for a new function-different from the function it had been selected for. An oftenquoted example of such functional shifts is the case of feathers: Initially, they might have evolved for temperature regulation, but, later, they were used for facilitating the flights (and water diving) of birds (50). In contrast, adaptation (from Latin adaptāre: ad-, toward, + aptāre, to fit) describes the emergence of a biological trait as a direct response to an environmental pressure. (54,55). Indeed, viruses, and, in particular, RNA viruses, seem well suited to test the predictions of simple evolutionary models: Their growth is fast and can be quantified, their mutation rate is high, and their population size can be easily controlled. It seems, therefore, that we could use the results of such laboratory experiments to at least qualitatively assess the validity of the predictions of our model. To do this, it is useful to consider the fluctuation relation derived in our restricted model.
We argued above that, under the assumptions of the restricted model, the probability of observing any population composition approximates, at long times, a Boltzmann-like distribution, Eq. 12. This distribution is ruled by the evolutionary potential Eq. 6 containing the φ n term. The fluctuation relation that ensues, Eq. 15, shows that the robustness to variations of the population increases as a function of the product μN (see also refs. 29 and 48). Most importantly, this relation connects the fluctuations of sensitivity, var{Ω n }, to the changes of the mean sensitivity as a function of the variation rate μ, d Ω n /dμ. The fluctuation relation further implies that, when the sensitivity to variations saturates and approaches a minimum, its fluctuations decrease. Numerical simulations qualitatively verify these predictions; see Box 2 and SI Appendix, section V.
Although we derived the enaptation of robustness and the fluctuation relation using the model [5], these results are expected to be approximately valid for any evolutionary stochastic dynamics appropriately including reproduction, variation, and selection [e.g., the aforementioned quasispecies as well as diffusion models (10,17)]. It is, indeed, these three elements that are essential here.
Possible Quantitative Experiments. Simple predictions of the fluctuation relation seem to be consistent with experimental results involving viral evolution. For instance, in ref. 56, two sets of viral populations founded from the same lineage were both evolved under strong selection but with different effective mutation rates μ. This was achieved by controlling the level of coinfection, which, in turn, through complementation, changed the effective probability of mutations μ. When subjected to mutation accumulation experiments, viruses evolved under low coinfection (i.e., under high μ) not only showed lower sensitivity to mutations, but also the fluctuations of sensitivity seemed to be significantly lower; see SI Appendix, section VIII. From the perspective of our fluctuation relation, this is a sign that the sensitivity has approached the saturation value. Although this qualitative agreement is encouraging, unfortunately, these experimental results do not present us with enough statistics to compare them quantitatively with predictions of evolutionary models. (e.g., the difference in the mean values of the growth rate between high-and low-mutation rate strains in ref. 56 is less than 1 SD of the variations in measurements). However, these interesting experimental results indicate that quantitative comparisons should be possible in principle. Other kinds of experiments which could provide quantitative data to be compared with the predictions of the RSV dynamics are 1) experiments on directed evolution of RNAs (57), proteins (58), or microbial communities (59); and 2) mutation accumulation experiments involving model organisms (11,12). If one is interested in the evolution of naturally replicating entities, however, viruses (including bacteriophages) are arguably the best candidates for relevant quantitative experiments.
Simplistic Nature of the Model. We conclude with a few additional comments about our model. It is clear that the dynamics described by Eq. 5 is a strong simplification of real evolutionary dynamics. Three major and related aspects have been here neglected: 1) For the sake of simplicity, we here disregarded direct interactions among organisms. De facto, organisms interact only indirectly through the environment, which affects their reproduction rate, variation probability, and overall population size. We refer to SI Appendix, section IB for the generalization of our model that accounts for interactions. 2) We deliberately focused on constant environments, while real environments do change and fluctuate, and their variations are affected by the population dynamics itself. Although such environmental changes can be formally incorporated in our model, we did not find them relevant at this stage of model development. Environmental changes could be described as an additional stochastic force: It would act alongside the forces described by Eq. 2, and it might favor types that quickly adapt to environmental changes (60)(61)(62). However, such additional forces depend on how the environment alters the population's reproduction, and hence are of a more idiosyncratic nature. Classification of "adaptation strategies" could form a complementary approach to studying the influence of changing environments (63,64). 3) We also disregarded any distinctions between genotype, epigenotype, and phenotype while introducing variations. This was done for the sake of generality. However, such distinctions definitely become important for changing environments. What differentiates genomes, epigenomes, and phenotypes (other than the molecular nature of their realizations) are the characteristic time scales at which they vary. These time scales reflect typical time scales over which information gathered about the environment remains useful. For constant environments considered here, the same information can be used forever, and there is no real need to differentiate between different time scales, and thus between different classes of variations.
In conclusion, we hope that, despite its simple nature, our approach may be of general interest. By connecting evolutionary dynamics to nonequilibrium statistical mechanics, it can render more precise some widely used notions, such as "evolutionary forces" or "chance" in evolution.
Data Availability. The code used for numerical simulations, as well as raw data are publicly available at GitHub, https://github.com/rjku/DarwinianEvolutionary Dynamics. All study data are included in the article and/or SI Appendix.